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Abstract 

Although Kolmogorov-Smirnov (KS) statistic is a widely used method, some weaknesses exist in investigating abrupt 
Change Point (CP) problems, e.g. it is time-consuming and invalid sometimes. To detect abrupt change from time series fast, 
a novel method is proposed based on Haar Wavelet (HW) and KS statistic (HWKS). First, the two Binary Search Trees (BSTs), 
termed TcA and TcD, are constructed by multi-level HW from a diagnosed time series; the framework of HWKS method is 
implemented by introducing a modified KS statistic and two search rules based on the two BSTs; and then fast CP detection 
is implemented by two HWKS-based algorithms. Second, the performance of HWKS is evaluated by simulated time series 
dataset. The simulations show that HWKS is faster, more sensitive and efficient than KS, HW, and T methods. Last, HWKS is 
applied to analyze the electrocardiogram (ECG) time series, the experiment results show that the proposed method can find 
abrupt change from ECG segment with maximal data fluctuation more quickly and efficiently, and it is very helpful to 
inspect and diagnose the different state of health from a patient's ECG signal. 
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Introduction 

Detecting abrupt change from time series, called CP detection, 
has attracted considerable attention in the fields of data mining 
and statistics. CP detection [1,2,3,4,5] has been widely studied in 
many real-world problems, such as atmospheric and financial 
analysis [1], intrusion detection in computer networks [2], signal 
segmentation in data stream [3], as well as fault detection in 
engineering systems [4,5]. A good method of CP detection is by 
comparing probability distributions of time series samples over 
past and present intervals [6,7], in which a typical strategy is to 
trigger an alarm for a CP as two distributions are becoming 
significandy different. Various methods of change detection follow 
this statistical framework, including the CUSUM (cumulative sum) 
[7], the GLR (generalized likelihood ratio) [8,9] and the change 
finder [4]. Generally, these approaches are limited by relying on 
pre-specified parametric models such as probability density 
models, autoregressive models, and state-space models. Therefore, 
these methods tend to be less flexible in real-world CP detection 
problems. 

In community of statistics, some non-parametric approaches for 
CP detection have been explored, in which non-parametric 
density estimation is used for calculating the likelihood ratio 
[10,11]. However, this kind of estimation is a hard problem 
[12,13], and may not be promising in practice. As a nonpara- 
metric method, the KS statistic quantifies a distance between the 
empirical distribution function of the sample and the cumulative 
distribution function of the reference distribution, or between the 
empirical distribution functions of two samples [14,15]. Moreover, 
the KS test is for the equality of continuous, one-dimensional 



probability distributions that can be used to compare a sample 
with a reference probability distribution, or to compare two 
samples. The null distribution of this statistic is calculated under 
the null hypothesis that the samples are drawn from the same 
distribution or that the sample is drawn from the reference 
distribution. The two-sample KS test is one of the most useful and 
general nonparametric methods for comparing two samples, as it 
is sensitive to differences in both location and shape of the 
empirical cumulative distribution functions of the two samples. 
Recendy, non-parametric KS statistic and its modified version are 
broadly investigated on several of application fields. For example, 
the use of the KS statistic for testing hypotheses regarding 
activation in blood oxygenation level-dependent functional MRI 
data [16]; modeling the cumulative distribution function of rub- 
induced AE signals and quantifying the goodness of fit with the KS 
statistic, to offer a suitable signal feature for diagnosis [17]; abrupt 
change point (CP) detection from electroencephalography signal 
(EEC) [18], and gene expression time series database [19]. 

On the other hand, Wavelet Transform (WT) is another 
promising approach for CP detection. In the past decade, WT 
approach has emerged as an important mathematical tool for 
analyzing time series [20,21,22,23,24]. It has found applications in 
anomaly detection, time series prediction, image processing, and 
noise reduction [20,24,25,26]. In particular, wavelets can repre- 
sent general functions at different scales and positions in a versatile 
and sophisticated manner, so the data distribution features can be 
easily extracted from different time or space scales [26,27]. The 
heart of wavelet analysis is Multi-Resolution Analysis (MRA), by 
which a signal can be decomposed into sub-signals of different size 
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resolution levels [7,27,28]. The properties of wavelets such as 
localization, orthogonality, multi-rate filtering, are essential for 
analysis of non-stationary and transient signals. WT can represent 
a general function in terms of simple, fixed building blocks at 
different scales and positions. These building blocks are generated 
from a single fixed function called mother wavelet by translation 
and dilation operations [29,30]. In addition, Haar Wavelet (HW), 
as a simpler WT, owns some attracting features including fast for 
implementation and ability to analyze the local feature. HW is a 
very useful to find discontinuities and high frequency changes in 
time series, and a potential candidate in modern electrical and 
computer engineering applications, such as signal and image 
compression, as well as abnormality detection from time series 
[7,28]. 

However, most of these methods above are time-consuming, 
and not scalable to large-scale datasets due to their time 
complexity. Moreover, some of them, e.g., KS, is occasionally 
insensitive even invalid for less significant data fluctuation, 
especially when abrupt change occurs near two endpoint areas. 
To detect abrupt change from time series quickly and 
efficiently, a novel non-parametric method is proposed based 
on multi-level HW and a modified KS statistic. In this method, 
we combine superiorities of both KS statistic and HW 
methods, and try to find an abrupt change in terms of 
maximal data fluctuation existing between two adjacent 
segments of a diagnosed time series. This paper is organized 
as follows. Section II implements the integrated HWKS 
method in detail. First, the two BSTs termed TcA and TcD, 
are constructed by means of multi-level HW from a diagnosed 
time series. Then, the framework of HWKS method is 
implemented by introducing a modified KS statistic and two 
search criteria based on TcA and TcD. Last, two HWKS- 
based algorithms are designed to implement CP detection from 
the diagnosed time series. Section III evaluates the perfor- 
mance of HWKS by comparing KS, HW, and T methods, via 
simulated time series and real ECG datasets. Section IV gives 
conclusion from previous sections. 

Method 

The flow diagram of the integrated HWKS framework 
(Fig. 1) is composed of three parts. First, the two BSTs, TcA 
and TcD are constructed from a diagnosed time series. 
Second, abrupt CP is detected from root to leaf nodes of 
TcA in terms of a modified KS statistic and two search rules. 
Last, the performance of HWKS is evaluated by comparing 
with KS, HW, and T methods. 

A. Construction of TcA and TcD 

Like all wavelet transforms, multi-level HW decomposes a 
discrete signal into two sub-signals with half its length. One sub- 
signal is a running average or trend; the other sub-signal is a 
running difference or fluctuation. HW is performed in several 
stages or levels [28] . It can be described using scalar products with 
scaling signals and wavelets. The discrete signals are synthesized 
by beginning with a very low-resolution signal, successively adding 
on details to create higher resolution versions, and ending with a 
complete synthesis of the signal at the finest resolution. Generally, 
by using A-level HW, a discrete time series signal Z = {z 1; z 2 ,..., 
Zfii}, can be decomposed into the £ th -level trend cA A , and k level 
fluctuations, i.e.,cD', cD 2 ,..., cD* k= 1, 2, .., log 2 N . As shown in 
Fig. 2, the Mevel HW is the mapping H k defined by [8], 



Z ^(cA^cD^cD^-'l.-lcD^cD 1 ), (1) 

The multi-resolution analysis (MRA) is the heart of 
wavelet analysis [23,29], in terms of MRA, we can conceptu- 
alize the process of HW as a projection of time series with 
size N to total jV different vectors zi; and w-„ termed as 
scaling signals and wavelet basis vectors, respectively. 
The discrete signal Z, average and detail signals are expressible 
as: 

k 

Z = A k + ^D',l<£<log 2 A, (2) 

N/2 k N/2 k 

A k = (Z- V k )V k = J2 (Z-v k )v k = £ (cA^v* (3) 

i=l 1=1 
N/2' 

D i = (ZW i )W i =^(cD it/ )w;, (4) 

7=1 

Thereafter, the following equations can be obtained, 
k 

Z = A k +J]D i 
1 = 1 

N/2 k k N/2' 

= £(z-v k )v k +x:£(z-w;)w; 

,=1 '=1 7=1 

N/2 k k N/2' 

= ]>>A A ,)v k +]r£(cD, i/ )w; 

1=1 i=l 7=1 

k 

= cA k -V k +^cD i -W J , 

7=1 

A k = cA k -V k 

= (cA k>I ,cA k>2 ,... > cA k>N/2 k)-(vf,v5,...,v^ /2k ) J (6) 

( n k _ k „k „k \ 

— V*l> d 2> • • • ' d N-l > a N>J 

D k = cD k W k 

= (cD M ,cD M ,...,cD kN/2k )-(w k ,w k ,...,w k /2k ), (7) 

= (d k ,d 2 k ,...,d k _ 1 ,d k ,) 

where v k is £-level Haar scaling signals, w k is A-level Haar 
wavelets, |v k | = |w k | =jV. 

In addition, we can represent HW with Mevel approximation 
and detail coefficient vectors by the following matrices, namely 
McA and McD: 
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based on TcA and modified KS statistic. 
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The performance evaluation of 
HWKS 
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Figure 1. The integrated scheme of HWKS method for fast CP detection, which includes three parts: construction of two BSTs, 
namely TcA and TcD; CP detection of HWKS in terms of two search criteria; and evaluation of HWKS method. 

doi:1 0.1 371 /journal.pone.0093365.g001 



McA = 



McD = 



cA 



i,i 



cA. 



cA 



l.N/2 

0 
0 

cD ljN/2 



cA mJ 0 

cD u .. 

... cD kJ 0 

cD m ,i 0 0 



(8) 



where 0<k<m = log 2 N, \<j<N/2 k . Suppose the size of a 
diagnosed sample Z is divisible k times by 2, we can further denote 
the / h element in cA k and corresponding averaged signal in A k , as 
well as the/ 11 element in cD k and corresponding detail signal in D k 
as follows: 



cA kJ = 



(V2)'k 



(9) 



cA* 



cA 3 



cA 1 



cA' r - 



c.V 



!cD l 



cD 

cD 3 
cD ; 
cD 1 



Figure 2. The diagram of multi-level HW for time-series signal 
Z, it is composed of Ar-level cA and cD vectors, i.e., the average 
and difference coefficients vectors. 
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= ak < 6= (V2) CAkj ' 



CD k j : 



dk,a = 



= dk ' b= (V2) CDk ' J ' 



(10) 



(11) 



(12) 



where 2<k<\og 2 N, and 2 k (J- 1) + 1 < i<j * 2 k ; a = 2 k (J-\) 
+ 1, c= 2 k (j — 1) + 2 (k ~ l \ and b = 2 k * j. Therefore, a diagnosed Z 
can be decomposed into cA and cD matrices by means of Mevel 
HW. Thereafter, as shown in Fig. 3, TcA and TcD are built in terms 
of McAand McD, as well as original elements in Z= {zi, Z-2,---, Zn}- 
In TcA, and TcD, non-leaf nodes in different level are constructed 
from McA, and McD, respectively; and then leaf nodes are derived 
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Figure 3. The diagrams of TcA and TcD, derived from a diagnosed time-series Zby means of /(-level HW. 

doi:10.1371/journal.pone.0093365.g003 



directly from the elements in Z. The constructions of TcA and TcD 
are implemented by Algorithm 1 in detail. 
Algorithm 1. 

Input: Z = {z ; : 1 < i < JV},a diagnosed time series 

Output: two BSTs, TcA and TcD 

Initiate the level of HW, k, (1 < k < log2(JV)); 

Declare two matrices, McA and McD; 

For i = 1 to £ do 

[cA ; , cD ; ] = Call HW(Z, t); 

McA(z) = cA ; ; McD( I ) = cD ; ; 
end 

Construct TcA and TcD from McA and McD, as well as Z. 
Output TcA and TcD 



D„ m {x)A{ ) 1/2 sup I G„(x) - F m (x) 
m + n .ver 



(14) 



if a change point c occurs in Z, there exists a value z (: satisfies 
F m (z c ) # G„(z c ), and D nm (z c )>5, i£[z,, z n ], (5ei?. 

As hypothesized -F„,(x) and G„(x) are not available, but instead, 
the e.c.d.f of F m (x) and G„{x) can be derived from two time series 
X and Z. Then, F m (x) and G„(x) are defined by, 



F m {x) = P m (X < x) = - V I(x, < x), 

111 ^ < 



(15) 



B. HWKS method based on a modified KS statistic 

KS test is one of the most useful and general non-parametric 
methods for comparing two samples, because it is sensitive to those 
differences in both location and shape of the empirical cumulative 
distribution functions (e.c.d.f) of two samples [31,32]. Suppose 
Y = {y\,...,ytf\ is a time series, we observe, 

Y=f(i/N) + X, i=l,...,N, 

where X = {x,} 1 = 1 N are discrete and centred i.i.d. random 
variables, and / is a noisy signal with unknown distribution. 
Thereafter, we can deal X as a normal time series with distribution 
function F m (x), and Tas an abnormal time series with distribution 
function G n (x). Then, we can assemble a diagnosed time series Z, 
and define it as below: 

Z = {X, Y} = {Z\,Z 2 } = {z\, . . . ,z c ,z c+ \, . . . ,z„}, (13) 

To detect an abrupt CP from Z, a modified KS statistic is defined 
to evaluate the distribution distance between Xand Z [15,33,34]: 



G n (x)=P n (Z<x)=~J2l(Zj<x), (16) 
7=1 

where F m (x) and G„(x) count the proportion of the sample points 
below level x. For any fixed point xeR, the law of large numbers 
implies that 

F m (x)=-J2l(x I <x)^E I( Xi <x) = F(x), (17) 



G„(x)= -^/(z,<x)^E I(z j <x) = G(x), (18) 

where F(x) and G(x) are true underlying distribution of two time 
series Xand Z> i-e. the proportion of the sample in the set (— 00 ,x] 
approximates the probability of this set. It is easy to say that this 
approximation holds uniformly over all xeR : sup \G n (x) — 

F m (x)\^0, the largest difference between F m (x) and G n (x) 
goes to 0 in probability. The key observation in KS test is that the 
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distribution of this supremum does not depend on the 'unknown' 
distribution P in diagnosed Z, if P is continuous distribution. In 
addition, we have, 



P((— ) 1/2 sup \G n (x)-F m (x)\ <S)- 
m+n xeR 

7/(c5)=l-2^(-l)'- 1 e 



-2i L a 



(19) 



where H(a) is the c.d.f. of KS distribution [35] . 

In terms of the definition of non-leaf nodes in TcA in equations 
(9), we can denote the/ h element Zj of Z and define a new element 
Zkj used in HWKS method as, 



z j = a Ki + E d 'J^ z J - ak J- 



(V2y(k+i) 



cA k 



2 k -Hy 

(V2Y(k+l) 



(20) 



(21) 



2 k 



-a k , a = --- = 



2 k 



ak,b<^av,Zb, 



where cA k j is a non-leaf node of TcA, z^j is a new element 
defined in terms of cA k j, and l</f<log 2 N, \<j<N /2 k ; 
a = 2 k (j— 1)+ 1, and b = 2 k * j. Then, a revised KS statistic for 
HWKS is defined as: 



D' mn (kj)M— —f' 2 sup \G n {z k j)-F m (z kj )\ 



nm 1/2 

— — ) sup 

n + m ZkjeR 



Y n | m 



7=1 



(22) 



where D' mn (kJ) measures the distribution distance between Xand 
Z at a selected node cA k j in TcA, and larger value of D' mn {kj) 
means that more significant change occurs in Z. 
Thereafter, we define a KS test for X and Z as, 



where c depends on the level of significance a, and can be found by 
using KS distribution when n is large. 



a = P(<5 # H 0 \H 0 ) = P(D' mn > c\H 0 ) « 1 - H(c), 



(27) 



If supsup|Z>' m „(/c,/)| < C(a), // 0 is true i.e., no change point 

k j 

occurs. On the other hand, if sup SWp\D' mn {kj)\ > C(a), then 

k j 

hypothesis H x is true i.e., an abrupt change is detected. 

C. HWKS-based CP Detection 

To detect abrupt change from diagnosed Z, an optimal path 
needs to be obtained from root to leaf nodes in TcA accurately and 
quickly. Therefore, as shown in Fig. 4, two search criteria are 
introduced in terms of TcA and TcD. The first search criterion is 
defined based on TcA as follows: 

Criterion 1. Suppose the current non-leaf node in TcA we 
selected is cAkj, and its left-child and right-child nodes are 
cA k _i 2j-i and cA k _i 2j, respectively, 

(a) if(D' mn (k- 1,2/- \)>D' m „(k- 1,2/) and D' m „{k- 1,2/- 1) 

> C(o<)) hold true, then the left-child node cA k _i 2j-i is 
selected to be involved into the current search path in TcA; 

(b) if (D' mn (k-l,2j-l)<D' mn (k-l,2j) and D' m (k- 1,2/) 

> C(o<)) hold true, then the right-child node cA k _i2j-i is 
selected to be involved into the current search path in TcA. 

Proof. In terms of the definitions of ztj and D' 
(21) and (22), D' m „(k- 1,2/- 1), and D' mn (k- 1,2/) can be written 

as 



D' m (k-l,2j-l)^^y< 2 \G n (.z k - laj - 1 )-F m (,z ll - 1 2 l - 



nm i /-> 

=<^) 1/2 

n + m 



n + m 



'\{^I{Za<zi-i,2j-\))-^'f^I(xb<Zk-\,y-\)} 

n a=\ m b=l 



(28) 



, (y/2)-«:- l) 
2* 



cA t . UH ))- ^E s 'W uj-,)} 



Hq ■ F m = G n vs. H 0 : F m ^ G n , 



(23) 



if null hypothesis is true then, the distribution of D' mn can be 
tabulated as it depends only on n. Moreover, if n is large enough 
then the distribution of D' mn is approximated by KS distribution. 
On the other hand, suppose F m ^ G n , since F m and G n axe the true 
c.d.f. of X and Z, bylaw of large numbers the e. c.d.f, F m will 
converge to G„, and for large n we will have, 



sup \G„(zkj)-F m (zkj)\>S,D' mn >(^) l ' 2 S, (24) 
z k,f R m+n 



If Hq fails then, 



D'mn>( ) 1/2 <5^ + O0 as W--- + 00. (25) 

m + n 



Therefore, to test Ho we make a detection rule, 



H, : D' 



• >c' 



(26) 



n+m 



«™ iraU x 5 - (V2)'(k-l) (V2TO-1) 

= (^,) 1/2 |{„(E^ * LJ 2 ^cA k . Uj)) - -g/fe < LJ^cA^)} 



where Z/t- 1,27-1 an d Z/fc-1,2/ ar e two diagnosed points in 



accordance with cA 



k-l,2j-l 



and cA k _! 2 j in TcA. In terms of 



Criterion 1, if (D' mn (k— 1,2/— \)>D' mn {k— 1,2/) holds true, as 
plotted in Fig. 5, it indicates that more significant distribution 
distance exists in the left sub-tree covered by cAv-i^j-i, than in the 
right one covered by cA k _i2j, and vice versa. That is, abrupt CP 
occurs in the left segment of Z with more probability than in the 
right one. On the other hand, If D' nm (k— 1,2/— 1)> C(a)) is 
satisfied, it means that the distribution distance overtakes a critical 
value given in an identical data distribution. 

Criterion 1 guarantees that if abrupt CP occurs in Z, the left or 
the right sub-tree with bigger distribution distance is selected to be 
involved into the current search path, and then the other half part 
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(a) Criterion 1 based on TcA 




Figure 4. The scheme of two search criteria for CP detection of HWKS. (a) Criterion 1 based on TcA, and (b) Criterion 2 based on TcD ensure 
that an optimal path of abrupt CP can be detected from root to leaf nodes in TcA. 
doi:1 0.1 371 /journal.pone.0093365.g004 



is discarded. Thereafter, an optimal search path can be detected 
from root to leaf nodes in TcA after about log 2 («) search steps. 
Unfortunately, if (D' m „(k — \ ,2j — \ ) = = D' mn (k — 1,2/) or 
(max(D'„ m (k — 1 ,2j — l),D' nm (k — 1 ,2jJ) < C(a)) are true, then 
Criterion 1 is invalid for CP detection. Therefore, another search 
criterion is introduced based on TcD, and defined as below. 

Criterion 2. Suppose D' mH (k— 1,2/— 1)= =H mn (k— 1,2/) 
or (ma.x(D' nm (k — 1,2/— \),D' nm {k— l,2/))< C(a)) is satisfied, 
the non-leaf node cD^j in TcD is selected, in accordance with 
the current non-leaf node cA^j in TcA, with its left-child node 
cDk— l,2j — l and right-child node cDj i _i i 2j, respectively, 

(a) if cD] c _i > 2j_i > cDt-i^j holds true, then the left-child 
node cAn^j-i is selected to be involved into the current 
search path in TcA; 

(b) if |cDk-i.2j-i | < |cDk-i,2j| holds true, then the right-child 
node cAk-i^j is selected to be involved into the current 
search path in TcA. 

Proof. In accordance with the definition of cD^j in equation 
(1 1), cD^.i 2j-l an d cD^-i^j can be written by 



cD k . UH = -( ZL "~ E zu) = DL(k-l,2j-l), (30) 

'" ■ ■ '-(, Lb = L x +\ 



I R l R 2 

cD k „ 1 , 2j = -( £ Z R"~ zZ ZRb) = D R (k-\,2jl (31) 

II Ra = R 0 Rb = R l +l 



where L 0 = 2 k (j- 1)+ 1, L, =2 k (j- l) + 2 (k - 2 \ L 2 =2 (k - l > 
(2/-1); R 0 =2 (k - 1 H2j-l) + l, J R 1 =2 (A - 1) (2/'-l) + 2 ( *- 2) , 
i?2 = 2^ k ~ v> (2j), and m = n = (y[Tf k In terms of equation 
(30), and (31), Z)l , and Dr can reflect data fluctuation of two 
segments, namely Z L and Z R in Z covered by cD k _i 2j-i and 
cD k _i2j, respectively. None loses of generalization, if 
|cDk_ i,2j-l | > | c Dk-i ,2j | is true, it means that bigger data 
fluctuation exists in Z.l than in Z.r. That is, a potential abrupt 
change exists in Z L covered by cD^i 2j-i with more probability 
than in Zr covered by cD^-i^j- Therefore, in terms of criterion 2, 
left-child node cA^.i^j-i is selected to be involved into current 
search path in TcA, and vice versa. 

In terms of two search criteria above, two HWKS-based 
algorithms are implemented to detect abrupt change from a 
diagnosed time series Z. In Algorithm 2, the distribution 
distance between X and Z is calculated in terms of a selected 
non-leaf node cAy in TcA. In this function, the confidence 
interval is set by a = 0.05, and C(cc)= 1.3258. For simplicity, we 
only output the nodes that the value of D' mn overtakes C(a), 
otherwise output zero. In Algorithm 3, an optimal path is 
detected from root to leaf nodes in TcA in terms of two search 
criteria, and then an estimated CP is obtained from a 
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Figure 5. The scheme of search strategy in Criterion 1 in terms of distribution distance between Xand Z. The dotted red line refers to a 
supposed CP, namely Z c , the solid black line, and green one stands for the points of z k . 12 j, and z k . 12j .i, respectively. 
doi:1 0.1 371 /journal.pone.0093365.g005 



Algorithm. 2. 

Input: a normal time series input X, and a diagnosed 
time series Z as well as current non-leaf node 
cA k j selected in TcA. 

Output: RS, the distribution distance between X and Z 
at the selected non-leaf node cA k j in TcA. 

Set Ca= 1.3258; 

RS = 0; 

Calculate Zkj ™ terms of cA k j; 

Sx = Calculate the e.c.d.f of z k j in X; 

Sz, = Calculate the e.c.d.f of z" kJ in Z; 

Sz 2 = Calculate the e.c.d.f of z k j in Z; 
Dl = abs(Sz r Sx); D2 = abs(Sz 2 -Sx); 
If (D1>D2) &&(D1>&) Then 
{Select Dl; 
RS = D1;} 
elseif (D 1 <D2) &&(D2>Ca) Then 
{Select D2; 
RS = D2;} 

elseif (Dl = =D2) | | (max (Dl, D2)<Ca) Then 

RS = 0; 
Endif 

Output RS(cA k j); 



Algorithm. 3. 

Input: X, Z TcA, and TcD derived from Z 

Output: The estimated abrupt CP from TcA and TcD. 

Set b = 1; JV= length (z); A = log 2 (JV); 

Set the first node of optimal search path is the root node in 
TcA: 

r]or i = 1 to k do 
a = k-i+l; 

Call Algorithm 2 to calculate the distribution 
distance between X and Z at tw° non-leaf node 
cA a _! 21,-1 and cA a _! 21,.! , respectively; 
Set SI =RS(cA a _i j 2b-i), S2 = RS(cA i ,i, 2b ); 
If (S1>S2) Then 
{The current selected node = cA a _i 2 b-ij 

b = 2b-\;\ 
elseif ((SKS2) Then 

{The current selected node = cA a _i 2 b', 
b = 2b;} 
elseif ((SI = =S2) Then 

{I%D a . li2b _i> cD a . li2b ) Then 

{ The current selected node = cA a _i 2 i,-i; 
b = 2b-l; } endif 
If(cD a .i :2 b-i< cD a _i i2 b) Then 



{ The current selected node : 
b = 2b; } endif} 
} Endif 
End for 

Output b, and Zb 



cA 



a- 1,2b) 



diagnosed Z. The pseudocodes can be found in Algorithm 2 
and 3 in detail. 

D. Evaluation of HWKS 

Many methods have been proposed for CP detection. In this 
part, the following typical methods are used to evaluate and verify 
the performance of the proposed HWKS method. 

KS statistic [15]. In KS method, firstly, we divide a 
diagnosed sample data Z into two segments, namely, Z m = {z b 
z 2 ,..., z m }, and zjv- m ={z m +i, z m+2 ,..., z/v- m }. Then, KS statistic 
for these two segments is defined by, 



D mn (x)A( 1 - ' ) 1 ' 2 sup \F N _ m {x)-FJx)\ 

xeR 

/at \ N-m m (32) 

= ( m(JV-m) )1/2sup| ^ m<x) _^ T{Zj<x% 

^ - veR , = ,„+! ,/=l 

where JV is the size of the diagnosed sample Z; m is the size of 
segments Zm> tna t is the current diagnosed position in Z. 

HW [28,33]. In this method, the fluctuation coefficient vector 
cD 1 ={cDij,cDi t 2, . ..,£■£>! W / 2 ,} is calculated by one-level HW 
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(A2) Z2 of size N=64, with k=9 
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Figure 6. The results of single simulation on single CP test position, with constant variance v=2, different sample size N, and CP 
test position k, by HWKS, KS, HW, and T, respectively. (A1 ) The estimated CP from diagnosed sample Z1 , with N = 32, k = 5. (A2) The estimated 
CP from diagnosed sample Z2, with N = 64, k = 9. (A3) The estimated CP from diagnosed sample Z3, with N= 128, k = 1 13. (A4) The estimated CP from 
diagnosed sample Z4, with W = 256, /c = 225. 
doi:10.1371/journal.pone.0093365.g006 



from a diagnosed time series Z, and then an estimated CP can be 



found by comparing the values of elements in vector cD 
elements in cD is defined as, 



/*2 



V2 



( £ (-D'" +1 z,), 



The 



(33) 



;=2(j-i)+i 



where 1 <j<N/2. To find an abrupt change, a critical value 5 is 
given in terms of an identical data distribution. If 

max IcD] ;| >\fl-& holds true, then an abrupt CP occurs in a 
\<j<N/T 1 

diagnosed sample Z. 

T-statistic[36]. T, also known as Welch'st-test, is used only 

when two population variances are assumed different (the two 

sample sizes may or may not be equal) and hence must be 

estimated separately. A diagnosed sample Z is divided into 

Z m ={z\, z 2 ,..., z m } and Z N . m = {z„ !+1 , z m+2 ,..., z N . m \. Then, T 

statistic is calculated as, 



Sf 



N — m 



(34) 



where Z m , and Zff _ m are the sample mean of two segments in 
diagnosed Z, S* is unbiased estimator of standard deviation, and m, 
n is the size of two segments in Z. 

Results and Discussion 

First, we evaluate the performance of HWKS on the simulated 
time series datasets, the sensitivity, efficiency, and accuracy of 
HWKS is analyzed by comparing KS, HW, and T methods. 
Then, we apply HWKS, and other three methods, to distinguish 
normal and abnormal ECG segments from the assembled ECG 
time series samples, and diagnose the different states of health 
from a patient's abnormal ECG time series segments. 



Table 1. The summarized results of single simulation with single CP test position. 





Z 


Size, N 


2 3 


2" 


2 5 


2 6 


2 7 


2" 


2 9 


2 ,o 




M 




CP, * 


2 


3 


5 


9 


113 


225 


449 


897 


Averaged 


HWKS 




Err 


0 


0 


1 


—3 


8 


-/ 


17 


-10 


5 






Acc 


1.0 


1.0 


.97 


.95 


.94 


.99 


.97 


.99 


.97 


KS 




Err 


-1 


-2 


5 


0 


6 


1 


-4 


-9 


3.5 






Acc 


.88 


.88 


.84 


1.0 


.95 


.99 


.99 


.99 


.94 


HW 




Err 


0 


1 


1 


17 


-51 


-51 


-93 


-503 


89 






Acc 


1.0 


.94 


.97 


.73 


.60 


.80 


.82 


.51 


.79 


T 




Err 


-1 


-2 


-1 


54 


-112 


-224 


62 


126 


72 






Acc 


.88 


.88 


.97 


.16 


.13 


.13 


.88 


.88 


.61 



doi:1 0.1 371 /journal.pone.0093365.t001 
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Figure 7. The results of multiple 600 simulations on single CP test position, with v= 2, different N, and k, by HWKS, KS, HW, and T, 
respectively. For samples with W= 32, k= 5; N = 64, k=9; N= 128, k= 1 13; and N = 256, /t = 225, (A1)-(A4) the distribution of e-CP, (B1)-(B4) the PDF 
of e-CP, and (C1)-(C4) the averaged e-CP, by HWKS, KS, HW, and T, respectively. 
doi:1 0.1 371 /journal.pone.0093365.g007 



Table 2. The summarized results of multiple 600 simulation on single CP test position. 







Z Size, N 


2 3 


2 4 


2 5 


2 6 


2 7 


2 8 


2 9 






M 


CP, k 


2 


3 


5 


9 


113 


225 


449 


897 


Averaged 


HWKS 


Tim 


.17 


.23 


.27 


.33 


.41 


.47 


.60 


.75 


.40 




Hit 


.39 


.15 


.10 


.06 


.03 


.03 


.01 


.01 


.10 




Err 


2 


2 


3 


3 


-3 


-3 


-5 


3 


3 




Acc 


.75 


.88 


.91 


.95 


.98 


.99 


.99 


.99 


.93 


KS 


Tim 


.09 


.15 


.36 


.68 


1.6 


3.4 


9.05 


23.8 


4.9 




Hit 


.0 


0 


.08 


.14 


.14 


.13 


.10 


.09 


.09 




Err 


-1 


-2 


-3 


-1 


-2 


-1 


-1 


1 


1.5 




Acc 


.88 


.88 


.91 


.98 


.98 


.99 


.99 


.99 


.95 


HW 


Tim 


.05 


.06 


.26 


.70 


1.15 


4.1 


8.7 


29.7 


5.6 




Hit 


0 


0 


0 


0 


0 


0 


0 


0 


0 




Err 


1 


4 


12 


25 


-46 


-95 


-186 


-380 


93.6 




Acc 


.88 


.75 


.63 


.61 


.64 


.63 


.64 


.63 


.67 


T 


Tim 


.6 


1.1 


2.2 


4.5 


9.2 


17.7 


36.7 


75.2 


18.4 




Hit 


.02 


.017 


.02 


.01 


.015 


.01 


.01 


.01 


.014 




Err 


0 


1 


2 


4 


-2 


-16 


-27 


-46 


12.3 




Acc 


1.0 


.94 


.94 


.94 


.98 


.94 


.95 


.96 


.95 



doi:10.1371/journal.pone.0093365.t002 
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Figure 8. The analysis of e-CP, hit rate, error and accuracy for multiple 1 00 simulations on different CP test positions, with different 
N, and k, by HWKS, KS, HW, and T, respectively. (A) The results of multiple samples Z1, with N= 16; (B) the results of multiple samples Z2, with 
N = 32; (C) the results of multiple samples Z3, with N = 64; (D) the results of multiple samples Z4, with N=128. 
doi:10.1371/journal.pone.0093365.g008 



A. Analysis on simulated time series 

In our simulations, the artificial time series are generated 
randomly in terms of normal distribution jV(0,l), i.e., jV[mean, 
M = 0, standard deviation, sd= 1), and then the normally distrib- 
uted datasets is used for abrupt change detection. Specifically, 
each diagnosed time series sample of size jVis composed of both a 
normal segment of size k, and an abnormal segment of size N-k, in 
which a designed abrupt change is contained by adding constant 
variation v to the normal random numbers of size N-k. 



Single simulation on single CP test position. First, CP 
detection is performed on single time series sample with single CP 
test position. The simulation results of CP detection are illustrated 
in Fig. 6, and the summarized analysis of error, and accuracy is 
shown in Table 1. For the proposed HWKS, it can detect the 
designed CP from samples of different sizes and CP test positions, 
with smaller error, and higher accuracy than HW, and T methods. 
For KS, it has the smallest averaged error in all four methods, but 
it has lower accuracy than HWKS, especially when sample size N 
is small. On the contrary, both of HW and T are worse, due to 



Table 3. The summary 


of mu 


Itiple simulations on different CP test positions. 










N. Z 
M \_ 


Size, N 


2 3 


2 4 


2 5 


2 s 


2 7 


2 s 


2 9 


2 ,o 


Averaged 


HWKS 


Tim 


.18 


.45 


.59 


.78 


.86 


1.18 


1.7! 


2.6 


1.04 




Hit 


.29 


.21 


.25 


.24 


.16 


.10 


.05 


.03 


.17 




Err 


1.25 


2.2 


3.2 


4.9 


6.8 


11.2 


21.6 


33.7 


10.6 




Acc 


.84 


.86 


.90 


.92 


.95 


.96 


.96 


.97 


.92 


KS 


Tim 


.14 


.50 


.91 


1.8 


4.1 


9.3 


22.9 


62.7 


12.79 




Hit 


.0 


.03 


.09 


.13 


.14 


.11 


.09 


.07 


.08 




Err 


3.5 


7.1 


8.8 


8.1 


9.5 


16.1 


31.0 


63.3 


18.4 




Acc 


.56 


.56 


.72 


.87 


.92 


.93 


.94 


.94 


.81 


HW 


Tim 


.03 


.18 


.45 


.77 


2.2 


8.2 


28.7 


72.7 


14.15 




Hit 


.08 


.08 


.04 


.02 


.02 


.03 


.02 


.03 


.04 




Err 


1.9 


3.3 


6.4 


13.2 


27.3 


61.8 


123.7 


251.9 


61.2 




Acc 


.76 


.79 


.79 


.79 


.78 


.76 


.76 


.75 


.77 


T 


Tim 


.78 


3.1 


6.1 


12.1 


24.3 


49.7 


101.7 


209.6 


50.9 




Hit 


.15 


.16 


.25 


.20 


.13 


.08 


.06 


.04 


.13 




Err 


1.38 


1.8 


2.1 


3.1 


5.7 


12.4 


25.3 


53.3 


13.1 




Acc 


.83 


.88 


.93 


.95 


.95 


.95 


.95 


.94 


.92 


doi:1 0.1 371 /journal.pone.0093365.t003 
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(A) Trend analysis for different sample size N=8 to 1024 



(B) Histogram analysis for the averaged time, hit rate, error and accuracy 
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Figure 9. The analysis of computation time, hit rate, error and accuracy on different sample size, for HWKS, KS, HW, and T, 
respectively. (A) The trend analysis for different sample size from W = 2 3 to 2 ,c , and (B) the histogram analysis for the averaged computation time, 
hit rate, error and accuracy. In (B), '1' stands for HWKS, '2' stands for KS, '3' stands for HW, and '4' stands for T. 
doi:1 0.1 371 /journal.pone.0093365.g009 



bigger error and lower accuracy than HWKS and KS. These 
simulation results indicate that HWKS has the best sensitivity and 
performance in four methods, especially HWKS is better than KS, 
due to smaller error and higher accuracy when abrupt change is 
located near the left or right boundary of samples with smaller size. 

Multiple simulations on single CP test position. Second, 
we test HWKS and other three methods by multiple 600 



simulations on single CP test position, with different N, and k. 
The representative results of CP detection are illustrated in Fig. 7, 
and the analysis of computation time, hit rate, error, and accuracy 
is summarized in Table 2. For the proposed HWKS, most of e- 
CPs are located near the designed CP position, with the shortest 
computation time, and the highest hit rate in all four methods, as 
well as smaller error, and higher accuracy than HW and T. For 
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Figure 10. The results of CP detection from assembled ECG time series of size N=2 k , k=9, 10, 14, by HWKS, KS, HW, and T, 
respectively. (A1)-(A6) the assembled ECG sample Z1-Z6; (B1)-(B6),(C1)-(C6) ,(D1)-(D6) ,(E1)-(E6) the e-CP detected from Z1 -26, by HWKS, KS, HW, 
and T, respectively; (F1)-(F6) the diagram analysis for the computation time, (G1)-(G6) the error of e-CP, and (H1)-(H6) the accuracy for Z1-Z6, 
respectively. In (F)-(H), 'V stands for HWKS, '2' stands for KS, '3' stands for HW, and '4' stands for T. 
doi:10.1371/journal.pone.0093365.g010 
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Table 4. The summary of CP detection from the assembled ECG samples. 



\. z 


Size, N 


2 9 


2 9 


2 ,o 


2 ,o 


2 11 


2" 


2 12 


2 12 


2 ,3 


2 ,3 


2 ,4 


2 14 




M 


CP, * 


50 


500 


300 


900 


600 


1400 


1600 


4000 


3000 


8100 


5000 


16000 


Averaged 


Time 


HWKS 


.032 


.036 


.032 


.033 


.035 


.034 


.034 


.034 


.041 


.036 


.060 


.049 


.038 




KS 


.014 


.017 


.037 


.038 


.116 


.118 


.395 


.397 


1.59 


1.46 


5.23 


5.58 


1.25 




HW 


.092 


.089 


.093 


.095 


.094 


.080 


.081 


.091 


.091 


.087 


.085 


.092 


.089 




T 


.074 


.079 


.149 


.141 


.296 


.293 


.718 


.696 


1.65 


1.66 


4.39 


4.38 


1.21 


Error 


HWKS 


7 


0 


0 


7 


7 


21 


0 


8 


8 


2 


12 


0 


6 




KS 


15 


227 


18 


15 


30 


28 


18 


28 


70 


6 


3 


15 


39.4 




HW 


2 


296 


96 


696 


396 


1196 


1396 


1744 


744 


3762 


662 


5092 


1340.1 




T 


49 


11 


299 


63 


599 


647 


2495 


95 


5191 


91 


11383 


83 


1750.5 


Accuracy 


HWKS 


.98 


1.0 


1.0 


.99 


.99 


.99 


1.0 


.99 


.99 


.99 
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KS, to some extent, it is better than HWKS and HW for higher 
accuracy and smaller error; however, it needs much more 
computation time, and has lower hit rate than HWKS when 
abrupt change occurs near the left or right boundaries of samples 
with smaller size JV. As for HW, it needs more computation time 
than HWKS and KS, and has the lowest hit rate and the biggest 
error in four methods. T is a method with the longest computation 
time in all four methods, and lower hit rate than HWKS and KS, 
although it has relatively higher accuracy than HWKS. These 
simulations show that, the proposed HWKS is a fast and efficient 
method, due to the shortest computation time and the highest hit 
rate in all four methods, as well as smaller error and higher 



accuracy than HW and T. In addition, HWKS has better 
sensitivity to less significant data fluctuation in sample with small 
size JV, than KS, HW, and T, especially when CP is located near 
the left or right boundary. 

Multiple simulations on different CP test 
positions. Third, for each diagnosed sample group, multiple 
100 simulations on different CP test positions are performed by the 
proposed HWKS and other three methods. In our simulations, for 
each JV, we select different 16 CP test positions from the different 
parts of diagnosed samples, i.e., the CP test position is designed by 
k = i * 2 (lo ^ N) -\ 7V>8, i= 1, 2, 3,..., 16. The selected results 
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Figure 11. The results of CP detection from abnormal ECG time series of size N=2 k , k="\0, 11, 15, by HWKS, KS, HW, and T, 
respectively. (A1)-(A6) the abnormal ECG sample Z1-Z6; (B1)-(B6), (C1)-(C6), (D1)-(D6), (E1)-(E6) the e.c.d.f derived from two segments ofZ1-Z6, by 
HWKS, KS, HW, and T, respectively; (F1)-(F6) the diagram analysis of the distance of e.c.d.f, and (G1)-(G6) the computation time of HWKS, KS, HW, and 
T in Z1-Z6, respectively. In (F)-(G), T stands for HWKS, '2' stands for KS, '3' stands for HW, and '4' stands for T. 
doi:10.1371/journal.pone.0093365.g011 
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Table 5. The summary of CP detection from abnormal ECG samples. 
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of simulations on single CP with different test positions are shown 
in Fig. 8, including e-CP, hit rate, error, and accuracy. In addition, 
the simulation results for different samples of size from 2' to 2 
are summarized in Table 3. In terms of computation time, hit rate, 
error, and accuracy, the trend analysis for different sample size JV, 
as well as the histogram of the averaged analysis for different 
methods are plotted in Fig. 9. 

For the proposed HWKS, it has the best performance for CP 
detection from samples Z of size JV from 2 3 to 2 7 , because of the 
shortest computation time, the highest hit rate, the smallest error, 
and the highest accuracy in all four methods. Meanwhile, HWKS 
is better than KS, HW, and T methods for samples Z of size JV 
from 2 to 2 , due to the shortest computation time, the smallest 
error, and the highest accuracy in all four methods, except that hit 
rate is slighdy lower than KS and T. For KS, in general, it has 
shorter computation time, bigger hit rate, smaller error and higher 
accuracy than HW. However, KS has the lowest hit rate, the 
biggest error, and the lowest accuracy as JV is below 2' in all four 
methods, which means that KS has worse sensitivity and 
performance for less significant data fluctuation, especially when 
JV is smaller. For HW, it generally takes shorter computation time 
than T, whereas, it has the smallest hit rate, the biggest error and 
lowest accuracy in all four methods. For T, it takes the longest 
computation time, although bigger hit rate, smaller error and 
higher accuracy than KS and T. These results show that HWKS 
has better performance and sensitivity to less significant data 
fluctuation near the left and right boundary of samples with 
smaller size JV and HWKS is an encouraging method for CP 
detection on simulated time series, due to shorter computation 
time, higher hit rate, smaller error, and higher accuracy than KS, 
HW, and T methods. 

B. Analysis on ECG time series 

To verify the performance of the proposed method further, we 
apply HWKS, and KS, HW, and T methods, to detect abrupt 
change from ECG time series provided by PhysioBank. In ECG 
experiments, we design the diagnosed ECG samples from different 
ECG datasets, including the MIT-BIH Normal Sinus Rhythm 
Database (NSRDB) [37], MIT-BIH Noise Stress Test Database 



(NSTDB) [38], and MIT-BIH Malignant Ventricular Arrhythmia 
Database (MVADB) [39,40]. 

CP detection from assembled ECG samples. First, we 
select a normal ECG dataset, 16265m from the NSRDB, and an 
abnormal ECG dataset, 118e()0m from the NSTDB, and then 
assemble the diagnosed ECG samples from different segments in 
the 16265m and 118e00m. Specifically, we take the normal ECG 
segment of size m as X m , and the abnormal segment of size n as 2~ n , 
respectively, and then assemble the diagnosed ECG sample 
Z= {X m , T n } — {xi,..., x m , yj,..., y n }. Meanwhile, we design 
another normal ECG segment from 16265m, i.e., X—{xi,..., 
x m+n}> as normal time series input. 

In this ECG experiment, a single CP test position is arranged 
near the left and right boundary of the assembled ECG sample. 
For the assembled ECG sample of size from JV=2 M to 2 1 * with 
different CP position k, the results of CP detection from Z1-Z6 are 
illustrated in Fig. 10, and then the analysis of computation time, 
error, and accuracy are summarized in Table 4. Comparing with 
KS, HW, and T methods, the results show that the proposed 
HWKS can estimate CP position more quickly, and distinguish the 
normal and abnormal segments from the assembled ECG samples 
more efficiently, with smaller error and higher accuracy. For KS, it 
has smaller error and higher accuracy than HW and T, whereas, it 
takes much more computation time than HWKS, HW, and T, 
especially when sample size JV gets bigger, meanwhile, KS is less 
sensitive for less significant statistic fluctuation, with bigger error 
near the right boundary. For HW, it is inefficient, because of the 
biggest error and lowest accuracy in all four methods. For T, it is 
also discouraging for longer computation time, bigger error and 
lower accuracy than HWKS and KS. Therefore, the proposed 
HWKS has the best performance in this ECG experiment out of 
all four methods. 

CP detection from abnormal ECG samples. To verify the 
performance of CP detection further, we apply the proposed 
HWKS, and KS, HW and T to analyze the abnormal ECG time 
series directly. In this part, we select the abnormal ECG segment 
from 118e00m in the NSTDB, i.e., Z= {yi,..., y n }, as a diagnosed 
ECG sample. Then, we take another normal ECG segment from 
16265m in the NSRDB, i.e., Z— {X n } = {xj,..., x n }, as normal 
input signal. To some extent, the distance of e.c.d.f can pardy 
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reflect the statistical fluctuation. Therefore, we take this variable as 
an indicator of the data fluctuation between two ECG segments 
divided by e-CP position. The results of CP detection, including 
the e-CP position, distance of e.c.d.f, and computation time, are 
plotted in Fig. 11, and summarized in Table 5. 

For abnormal ECG samples Z1-Z6 with different size jVfrom 
2 10 to 2 15 , HWKS can detect abrupt change position, and then 
divide the original ECG sample into two parts, with the shortest 
computation time out of four methods, and bigger distance of 
e.c.d.f than HW and T methods. For KS, it can detect CP with the 
maximal distance of e.c.d.f; however, it needs longer computation 
time than HWKS, HW and T, especially for ECG sample with big 
size JV. On the other hand, for HW, it is inefficient, due to smaller 
distance of e.c.d.f than HWKS and KS. T is also inefficient, 
because of longer computation time than HWKS, HW, the 
smallest distance of e.c.d.f, and the invalid e-CP position for most 
of the abnormal ECG samples. 

Especially, the results of CP detection from Z1-Z3 plausibly 
indicate that, a patient seems recovering from abnormal state of 
health, after overtaking the critical e-CP position detected by 
HWKS. On the contrary, the results from Z4-Z6 suggest that, a 
patient is encountering a risky situation from the former state of 
health, after going through the vital e-CP position. These results 
indicate that HWKS can capture abrupt change position from a 
diagnosed ECG sample quickly and efficiently, and the detected 
CP is very useful to find a critical time from ECG time series, 
where a patient might encounter an important conversion between 
two different states of health. Therefore, HWKS is an efficient and 
encouraging method for detecting abrupt change from abnormal 
ECG time series, and it is very meaningful in inspecting and 
diagnosing different states of health from diagnosed ECG time 
series more quickly and efficiently. 

Conclusion 

In this paper, based on HW and a modified KS statistic, a novel 
HWKS method is proposed for CP detection from large-scale time 
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